Following a comment by Meromorphic on the discourse here, this PR implements the matrix exponential function expm. It does so in a new stdlib_linalg_matrix_functions submodule which might eventually accomodate later on implementations of other matrix functions (e.g. logm, sqrtm, cosm, etc).

Proposed interfaces

• E = expm(A [, order, err])
• call matrix_exp(A [, order, err]) (in-place)
• call matrix_exp(A, E [, order, err]) (out-of-place)

where A is the input n x n matrix, order (optional) is the order of the Pade approximation, and err of type linalg_state_type for error handling.

Key facts

The implementation makes use of a standard "scaling and squaring" approach to compute the Pade approximation with a fixed order. It is adapted from the implementation by John Burkardt.

Progress

• Base implementation.
• Tests
• in-code Documentation
• Specifications
• Example

Possible improvements

The implementation is fully working and validated. Computational performances and/or robustness could potentially be improved eventually by considering:

• Automatic estimation of the best order for the Pade approximation based on the data. This is used for instance in scipy.linalg.expm.
• Balancing is not used currently, but it might improve the robustness.

In practice, it has however never failed me so far.

cc @perazz, @jvdp1, @jalvesz

PS : This is actually a re-opening of #968 which was automatically closed when I deleted my own fork of stdlib (don't ask me why, I can't remember).